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SECTION  I 
INTRODUCTION 

The  motion  of  a  rigid  body  about  a  fixed  axis  is  the  next  problem 
after  the  dynamics  of  particles  in  classical  mechanics.  The  rigid  body 
problem  contains  two  parts,  the  first  of  which  is  the  motion  of  the 
center  of  mass  of  the  body,  and  the  second  of  which  is  the  motion  of  the 
body  with  respect  to  the  center  of  mass.  It  is  the  second  part  which  has 
remained  an  unresolved  problem  in  classical  dynamics.  The  problem  can  be 
completely  specified  mathematically,  but  in  most  cases  the  solution  cannot 
be  expressed  in  closed  form. 

The  versatility  of  numerical  methods  along  with  the  availability  of 
digital  computers  has  revolutionized  the  study  of  many  problems.  Indeed, 
the  use  of  finite  element  methods  has  allowed  the  simulation  of  many 
difficult  nonlinear  problems  in  mechanics. 

The  following  subsections  contain  a  description  of  the  general  rigid 
body  problem  to  be  addressed  in  this  research,  as  well  as  a  brief 
historical  section  on  the  development  of  solutions  to  the  problem.  The 
objective  and  scope  of  the  research  is  stated  in  subsection  3.  The 
results  of  a  literature  search  are  found  in  subsection  4. 

1.  PROBLEM  DESCRIPTION 

Consider  a  rigid  body  of  arbitrary  shape,  moving  about  a  fixed  point 
within  the  body  due  to  the  effect  of  applied  torques  or  forces  (see 
Figure  1).  The  angular  momentum,  L,  of  the  body  is  given  by 

the  expression 

-  “  (1) 

L  =  Ito 

where  I  =  moment  of  inertia  of  the  rigid  body 
u  =  angular  velocity  of  the  rigid  body 
and  I  and  <o  are  for  convenience  referred  to  the  principal  axes  in  the 
body. 
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In  an  inertial  frame,  Newton's  second  law  for  rotational  motion 
is  written  as  follows 


=  N 


(2) 


where  N  =  vector  of  applied  torques. 

Equation  2  can  be  re-written  for  moving  coordinates  in  the  following  form 


(dLl 

[dLl 

[dtJ 

dt 

space 

V 

+  a)  x  L 


body 


(3) 


-Substituting  Equation  1  and  Equation  2  into  Equation  3,  one  arrives  at 
N  =  Iu>  +  u>  X  (Ill)) 

Equation  4  can  be  expanded  for  the  three  body  axes 


N1  = 

Vl ' 

u2- 

(5) 

n2  = 

1 2^2  * 

(V 

V“l“3 

(6) 

n3  = 

T3“3  ' 

(h  - 

1 2 ) 

(7) 

Equations  5-7  are  known  as  Euler's  equations  of  motion  for  a  rigid 
body  rotating  about  a  fixed  point.  The  solutions  w-| ,  w-j  describe 
the  motion  of  the  rigid  body  in  a  set  of  coordinates  fixed  in  or 
parallel  to  the  body. 

2.  BACKGROUND 

The  rigid  body  problem  of  subsection  1  has  well  defined  solutions 
when  the  body  is  symmetric,  i.e.,  two  or  three  moments  of  inertia  are 
equal.  However,  when  the  rigid  body  is  of  arbitrary  shape,  i.e., 
unsymmetric,  the  solutions  are  either  difficult  to  produce  or  unknown. 
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Classical  dynamics  texts  such  as  Goldstein  (Reference  1),  MacMillan, 
(Reference  2),  Whittaker  (Reference  3),  and  Sommerfeld  (Reference  4), 
give  detailed  accounts  of  special  cases  where  solutions  can  be  found. 
However,  these  solutions  are  in  the  form  of  elliptic  integrals  and  other 
special  functions,  and  require  clever  substitutions  in  order  to  manipulate 
the  problem  into  a  tractable  one. 

A  concise  summary  of  the  historical  development  of  solutions  to 
the  general  rigid  body  problem  can  be  found  in  Leimanis  (Reference  5). 

The  following  names  highlight  the  solution  of  the  general  rigid  body 
problem. 

1.  Euler  (1758):  Symmetric  body  with  no  applied  forces. 

2.  Lagrange  (1788):  Two  moments  of  inertia  equal  and  no  applied 
forces . 

3.  Kovalevskaya  (1888):  Two  moments  of  inertia  equal,  third 
moment  of  inertia  equal  to  half  the  other  two. 

4.  Staude  (20th  century):  llnsymmetric  top  under  force  of  gravity. 

5.  Bottema:  Stability  of  Staude  problem. 

6.  Grammel  (1948):  Use  of  approximating  functions  in  a  recursive 
scheme. 

A  brief  description  of  Grammel 's  work  can  be  found  in  Section  III. 

3.  OBJECTIVE  AND  SCOPE 

The  problem  outlined  in  subsection  1,  as  well  as  the  known  solutions 
in  subsection  2,  indicates  the  need  for  an  approximation  method  based  on 
modern  digital  computer  methods.  The  need  to  study  the  motion  of 
unsymmetric,  or  nearly  symmetric  rigid  bodies  might  result  from  the  design 
of  gyroscopes,  an  unsymmetric  mechanism,  or  a  satellite  antenna.  Thus, 
the  objective  of  this  research  is  to  investigate  the  use  of  an  approxi¬ 
mation  method  to  obtain  the  description  of  motion  of  a  general, 
unsymmetric  rigid  body  under  the  Influence  of  general  torques  or  forces. 
The  prevalence  of  finite  element  methods  in  other  areas  of  mechanics 
provides  a  motivation  to  use  a  trial  function  technique  analogous  to 
finite  element  methods  to  approximate  the  motion  of  a  rigid  body. 
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A  secondary  objective  is  to  define  a  method  which  is  inherently  easy  to 
use  and  not  highly  dependent  on  special  knowledge  or  manipulations. 

The  scope  of  this  research  is  limited  to  a  comparison  of  weighted 
residual  methods  (discussed  in  Sections  III  and  IV)  with  several  well 
known  numerical  methods.  Thus  the  purpose  of  the  work  is  concerned 
with  computer  approximation  methodology  to  generate  numerical  solutions 
to  the  general  rigid  body  problem  of  subsection  1. 

4.  LITERATURE  SEARCH 

A  1 iterature. search  was  conducted  in  order  to  identify  any  existing 
work  in  the  area  of  approximate  solutions  in  rigid  body  rotational 
dynamics.  The  University  of  Dayton  file  of  theses  was  inspected,  and  no 
work  had  been  done  on  this  problem.  A  search  was  made  of  books  on 
finite  element  methods.  These  sources  included  the  references  by 
Prenter  (Reference  6),  and  Zienkiewicz  (Reference  7),  as  well  as  the 
SAE  proceedings  (Reference  8).  The  material  addressed  under  rigid  body 
dynamics  in  finite  element  texts  was  generally  concerned  with  the 
determination  of  the  vibrational  modes  of  rigid  links,  bars,  or  other 
members,  and  as  such  was  not  relevant  to  the  issues  in  this  research. 

A  computerized  search  was  done  by  way  of  the  Lockheed  data  system  at  the 
Air  Force  Wright  Aeronautical  Laboratories  (AFWAL)  library,  wherein  the 
following  files  were  examined: 

Compendex,  Engineering  Index 

Dissertation  Abstracts 

ISMEC  (Info.  Service  in  Mechanical  Engineering) 

NTIS  (National  Technical  Information  Service) 

DDC  (Defense  Documentation  Service,  government  sponsored  research) 

NASA 

The  results  of  the  above  searches  were  a  series  of  computer  printouts  with 
finds,  each  containing  titles,  authors,  locations,  and  abstracts.  There 
were  a  few  listings  which  led  to  further  investigation,  which  are  listed 
in  the  List  of  References  and  the  Bibliography.  For  example,  the  papers 
by  Likins  (References  9-12),  Kraige  and  Skaar  (Reference  13),  and 
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Huston  and  Passerello  (References  14-16),  were  read  for  relevance  as  a 
result  of  the  searches.  In  general,  the  listings  referred  to  either 
multiple  rigid  body  analysis  (see  Section  III),  or  to  finite  element 
analysis  of  rigid  plates,  shells,  or  beams.  No  source  was  Identified 
which  duplicated  the  scope  of  the  project. 
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SECTION  II 

A  REVIEW  OF  METHODS  OF  DYNAMICS  AND 
SOLVABILITY  CONSIDERATIONS 


The  purpose  of  this  chapter  is  to  briefly  review  three  widely  used 
methods  in  dynamics  and  to  show  their  similarities.  These  include 
the  methods  of  Euler,  Lagrange,  and  Hamilton.  From  these  the  problem 
in  subsection  1  of  Section  I,  is  described  and  a  choice  of  equations  is 
made.  In  subsection  5  of  this  section  a  solvability  consideration  is 
described  such  that  bounded,  reasonable  approximations  can  be  expected. 


Eulerian  dynamics  is  concerned  with  the  conservation  of  angular 
momentum,  and  can  be  thought  of  as  another  statement  of  Newton's  second 
law  for  rotational  dynamics.  Lagrange's  equations  form  an  elegant 
alternative  to  the  Euler-Newton  method,  making  use  of  kinetic  and 
potential  energy  statements  and  a  function  of  generalized  coordinates 
and  velocities.  In  Hamilton's  method,  the  n  Lagrange,  second-order, 
differential  equations  are  replaced  by  2n  first-order,  partial  differ¬ 
ential  equations. 

1.  EULER'S  EQUATIONS 

Euler's  equations  for  the  rigid  body  were  given  in  Section  I  by 
Equations  5-7,  i.e., 

N1  =  Vl  "  ^2  " 


“  I 2^2  —  ( 


N3  =  "  ^1  " 


These  equations  specify  the  angular  velocity  of  the  rigid  body  about  its 
center  of  mass.  The  variables  are  in  body  coordinates.  The 
equations  are  coupled  and  nonlinear.  They  are  coupled  in  that  each 

variable  u>.  occurs  in  equations  other  than  the  one  containing  its 

*  ^ 

derivative  The  equations  are  nonlinear  in  that  each  equation  contains 
a  product  term  The  solution  of  the  equations  is  complicated  by  this 
nonlinear  coupling. 
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2.  LAGRANGE'S  EQUATIONS 

In  the  formulation  of  Lagrange's  equations,  the  crucial  problems  are 
finding  the  kinetic  and  potential  energy  expressions,  solving  the  dif¬ 
ferential  equations,  and  handling  of  any  nonholonomlc  constraint  forces. 
Lagrange's  equations  are  given  by 

d  3L  3L  _  nj 

™  94j  ~  3qj  (8) 

where  L  *  Lagranglan  *  T-V 
T  *  kinetic  energy 
V  =  potential  energy 
q.  =  generalized  coordinates 

J 

cjj  =  generalized  velocities 
■  generalized  forces 

The  generalized  coordinates  q^  do  not  have  to  be  the  specific  space 
coordinates  of  the  problem,  but  can  be  "quasi -coordinates" ,  l.e., 
coordinates  suited  to  the  problem  at  hand  and  not  of  necessity 
Integrable  combinations  of  the  time  derivatives  of  other  coordinates. 

The  quasi-coordinate  concept  is  explained  in  Meirovitch  (Reference  17, 
pages  137  and  157-160). 

Lagrange's  equations  for  quasi -coordinates  for  the  rigid  body 
problem  are  given  by 

3f  l£l  +  !“1  [£l  *  tN]  (9) 

where  [w]  is  the  vector  of  the  coordinates  ,  o>2,  and  [N]  is  the 
vector  of  generalized  torques.  The  angular  velocities  ,  w2,  0)3  are 
not  Integrable  time  rates  of  change  of  angular  displacements,  but  of  the 
quasi -coordinates.  The  advantage  In  using  this  formulation  Is  that  the 
equations  are  written  In  terms  of  an  orthogonal  set  of  axes. 
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3.  HAMILTON'S  EQUATIONS 

Hamilton's  equations  are  more  general  than  the  methods  of  Euler  and 
Lagrange.  The  Hamiltonian,  H,  Is  related  to  the  Lagranglan,  L,  by  the 
relationship 

H  =  Pa  q“  -  L  (10) 

where  p^  =  generalized  momentum  components 
qa  =  generalized  velocities 
Hamilton's  equations  are  given  by 


daa  = 

3H 

(ID 

dt 

3p 

a 

d£a  = 

-3H 

(12) 

dt 

3qa 

The  H  is  expressed  as  functions  of  the  coordinates  q“  and  momenta  p“, 
whereas  in  Lagrangian  dynamics  L  is  expressed  in  functions  of  the 
coordinates  and  velocities.  Hamilton's  equations  for  the  rigid  body 
problem  can  be  found  in  Pars  (Reference  18)  and  Webster  (Reference  19). 

The  resulting  expressions  for  Equations  5-7  are  six  partial  differential 
euqations  involving  the  Euler  angles. 

Hamilton's  equations  do  not  facilitate  the  solution  of  particular 
problems,  but  lead  to  important  theoretical  generalizations  in  fields  such 
as  quantum,  statistical,  and  celestial  mechanics.  There  is  no  general 
technique  for  solution  of  the  equations  in  closed  form.  However,  the 
resulting  first  order  differential  equations  are  sometimes  amenable  to 
solution  by  well  established  methods  in  specific  cases.  In  order  to 
handle  effectively  a  given  problem,  a  knowledge  of  canonical  or  contact 
transformations,  generating  functions,  and  Jacobi's  theory  is  required. 

4.  CHOICE  OF  METHOD 

It  can  be  shown  that  Euler's  equations  are  equivalent  to  Lagrange's 
equations  for  quasi-coordinates  (see  Melrovltch,  Reference  17,  pages 
157-169).  The  quasi-coordinate  approach  provides  the  most  direct  means 
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for  representing  a  general  rigid  body  of  arbitrary  shape  In  a  set  of 
orthogonal  body  coordinates.  Hamilton's  equations  were  judged  to  be  too 
complicated  to  be  appropriate  and  were  not  pursued  any  further. 


For  convenience.  Equations  5-7  are  re-written  Into  a  normalized 
form  as  follows: 

0>1  -  A“2t°3  = 

(*>2  ■  Bujto-j  =  M2 
0)3  -  Co)^u>2  =  M3 


(13) 

(14) 

(15) 


where 


M, 


M, 


M, 


V1!* 

A  3  (I2  - 

I2)/I 

n2/i2. 

B  3  (I3  - 

Ix)/I 

N3/!3, 

C  =  (Ij  - 

i2)n 

5.  SOLVABILITY  OF  THE  EULER  EQUATIONS 

In  light  of  the  non-linearity  In  Equations  13-15,  it  is  fundamental 
to  ask  whether  a  solution  exists  and,  if  so,  if  it  is  unique.  These 
equations  are  of  first  order,  and  there  is  much  theory  applicable  to 
first  order  systems.  See  for  example  Birchoff  and  Rota  (Reference  20, 
Chapter  6). 

Each  of  the  Euler  equations  can  be  re-written  In  the  form 


*  F(uj,  u>2»  Wj,  t)  3  M^  ^  ^2^^ 

u)2  3  G(o)j,  u>2»  u>2*  t)  3  M2  +  BoijUg 

0)3  ■  H(a>^,  ojg ,  0)3,  t)  ■  M3  +  Co> 

lu2 


(16) 

(17) 

(18) 
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The  critical  questions  consist  of  the  continuity  of  F,  G,  H  and  their 
partial  derivatives,  and  the  boundedness  of  F,  G,  H.  If  these  conditions 
can  be  inferred,  at  least  within  reasonable  intervals  of  computation, 
then  the  system  has  a  unique  solution.  One  means  to  test  for  boundedness 
Is  to  Introduce  the  Lipschltz  condition,  which  states  that  for  two  points 
x  and  y  on  a  region  R: 

|f(x,  t)  -  f(y,  t)|  <_  L | x  -  y| ,  (x,  t),  (y,  t)  e  R  (19) 

L  is  the  Lipschitz  constant,  which  Is  an  arbitrary  scalar.  Thus,  if  f  is 
bounded,  an  L  can  be  found  which  satisfies  a  Lipschitz  condition  on  an 
interval.  First  order  differential  equation  theory  reveals  that  if 

=  F((o-j ,  u2,  u>3,  t)  satisfies  a  Lipschitz  condition  on  the  domain  tj ,  t2 
then  there  is  at  most  one  solution  w^(t)  for  a  given  initial  condition, 
and  thus  uniqueness  can  be  inferred.  In  the  computer  simulations  in 
Section  IV,  this  property  is  used  to  monitor  the  bound  on  consecutive 
approximations.  An  unbounded  or  divergent  result  is  used  to  terminate 
the  problem. 
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SECTION  III 

REVIEW  OF  COMPUTATIONAL  METHODS 

The  objective  of  the  project  is  to  investigate  the  use  of  a  computer 
approximation  method  for  study  of  the  general  rigid  body  problem  of 
Section  I.  In  Section  II,  the  choice  of  equations  is  made.  The  desir¬ 
ability  of  a  trial  function  approach  Is  a  sub-objective,  and  because  of 
it  a  number  of  finite  element  methods  have  been  surveyed  for  relevance. 

The  weighted  residual  method  results  as  a  trial  function  approach. 

In  studying  the  applicability  of  weighted  residual  methods,  the 
need  is  also  generated  to  compare  the  results  with  another  known  method. 
Subsection  1  contains  a  description  of  weighted  residual  techniques.  In 
subsection  2,  several  known  numerical  methods  are  described.  Subsection  3 
provides  a  rationale  for  the  choice  of  methods  used  in  the  project.  In 
subsection  4,  a  brief  review  is  provided  of  Grammel's  trial  function 
approach,  as  well  as  multi-rigid  body  analysis. 

1 .  WEIGHTED  RESIDUAL  METHODS 

A  commonly  used  approach  in  approximation  of  ordinary  differential 
equations  is  the  weighted  residual  method.  This  method  is  characterized 
by  the  use  of  trial  solutions  with  undetermined  parameters,  which  are 
used  to  force  a  "residual"  to  approach  zero.  Four  common  weighted 
residual  techniques  are:  collocation,  subdomain,  Galerkin,  and  least 
squares.  The  application  of  these  techniques  is  illustrated  in  a  simple 
and  straightforward  manner  in  Crandall  (Reference  21). 

In  the  weighted  residual  method  a  trial  family  of  approximate 
solutions  is  selected.  For  example,  consider  the  problem 

x  =  f(x,  t),  x(0)  =  1  (20) 

One  could  select  the  polynomial  family 

x  =  1  +  Cjt  +  Cgt^  (21 ) 
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as  a  trial  function  with  undetermined  parameters  c^ ,  and  Cg.  The  trial 
function  is  chosen  to  satisfy  the  initial  conditions  of  the  problem 
independently  of  the  parameters  c-j  and  Cg.  The  residual,  R(t)  is 
formed  by  collecting  all  terms  of  the  governing  differential  equation 
together.  In  the  example, 

R  =  x  -  f(x,  t)  (22) 

If  the  trial  solution  is  the  exact  solution,  then  R  =  0.  The  four  tech¬ 
niques  are  different  criteria  for  determining  the  undetermined  parameters 
(Cp  C2  in  the  example)  over  each  interval. 

a.  Collocation 

With  the  collocation  method,  distinct  locations  throughout  the 
Interval  are  chosen  at  which  the  residual  is  set  to  zero.  The  number 
of  points  corresponds  to  the  number  of  undetermined  parameters  (two  in 
this  case).  The  unknown  parameters  are  found  by  solving  a  set  of 
simultaneous  equations  generated  from  the  values  of  the  residual  at 
the  specified  locations  in  the  interval.  The  choice  of  locations 
provides  an  effective  weighting  factor  on  the  residual.  It  is  possible 
to  design  an  algorithm  to  examine  and  adjust  these  weights  for  an  optimum 
between  residual,  computer  time,  etc.  The  interval  size  and  trial 
functions  can  also  be  adjusted  according  to  the  extent  of  non-linearity 
in  the  problem.  These  characteristics  can  be  used  in  the  remaining 
three  methods. 

b.  Subdomain 

In  this  method  the  interval  Is  divided  into  as  many  subdomains 
as  there  are  parameters.  The  Integral  of  the  residual  over  each  of  the 
Intervals  is  set  to  zero  to  provide  equations  for  determining  the 
unknown  parameters.  For  example,  using  the  subdomain  (0,  1/2)  and  (1/2,  1), 

*5 

/  Rdt  =  ' 

0 

1 

/  Rdt  =  0 

*5 


(23) 
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c.  Galerkin 

In  this  case,  weighted  averages  of  the  residual  over  the  entire 

interval  are  required  to  vanish.  The  weighting  functions  chosen  are  the 

2 

same  functions  of  t  used  In  the  trial  family  (i.e.,  t  and  t  ,  in  the 
example).  For  example, 

1 

/  tRdt  =  0 
0 

!  (24) 

I  t2Rdt  =  0 
0 


d.  Least  Squares 

In  this  case  the  criterion  Is  to  minimize  the  integral  of  the 
square  of  the  residual  over  the  Interval.  The  formulas  are  given  by 


f  R2dt  -/  R  ~  dt  -  0 
3cl  1  0  3cl 


h~  I  R2dt 
3c2  0 


1  3R 

/  R  ~  dt  •  0 
0  3C2 


(25) 


e.  General  Remarks 

The  four  methods  above  can  be  summarized  by  the  statement 

1 

/  WRdt  =  0  (26) 

0 

where  W  represents  a  weighting  factor.  The  choice  of  weighting  factors 
for  the  four  methods  discussed  are:  In  collocation,  the  Dirac  delta 
function  6(t);  in  subdomain,  the  unit  step  function;  in  Galerkin,  the 
trial  functions;  and  in  least  squares,  the  factors  fjr  ,  |jf  .  The 
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collocation  method  has  an  advantage  In  using  the  sifting  property  of  the 
delta  function  in  that  the  residual  can  be  evaluated  directly  without 
integration,  i.e., 

oo 

/  R ( t ) 6 ( t  -  tx)dt  =  R(tx)  (27) 

-oo 

Crandall  shows  an  example  problem  with  a  known  solution  and  com¬ 
pares  the  methods  with  a  Taylor  series  approximation.  In  that  case,  all 
four  give  significantly  better  results  than  the  Taylor  series,  with 
the  least  squares  and  subdomain  methods  the  best  over  the  interval 
chosen  in  the  example.  These  results  will  vary  depending  on  the  problem, 
the  interval  size,  and  the  trial  functions. 

2.  NUMERICAL  METHODS 

The  purpose  of  this  section  is  to  provide  a  brief  summary  of  the 
numerical  methods  which  are  considered  in  this  project. 

a.  Improved  Euler  Method 

The  basic  Euler  method  consists  of  the  construction  of  a  tangent 
line  approximation  on  a  point  by  point  basis  along  the  curve  f(x,  y). 

If  y  =  f(x,  y),  with  y(xQ)  =  yQ  then  the  Euler  method  consists  of 

Vl  =  yn  +  hf(x’  y)  =  yn  +  hyn  (28) 

where  n  =0,  1 ,  2,  ... 
and  h  =  step  size. 

The  Euler  method  is  very  easy  to  use;  however,  the  local  formula  error 

2 

is  proportional  to  h  ,  and  thus  the  method  is  not  very  accurate. 

In  the  improved  Euler  method,  the  above  formula  is  complemented  by 
the  average  of  the  values  at  two  points,  thus  improving  the  tangent  line 
estimate.  The  improved  Euler  formula  is  given  by 

yn  +  f[x  +  h,  y  +  hy  ]  h 

yn+l  =  yn  +  ~ - [L-1— 1 3 - -  (29) 

3 

The  error  Is  proportional  to  h  .  This  is  the  simplest  "predlct-correct" 
method. 
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b.  Three  Term  Taylor  Series 

The  Euler  formula  (Equation  28)  is  in  essence  a  two  term  Taylor 
series  approximation.  This  can  be  made  more  accurate  by  using  the  first 
three  terms.  Thus, 


Vl  *  +  <30) 

where  yR  =  f(x,  y) 

and  yn  =  fx(V  yn>  +  V  V  yn)yn 

The  local  formula  error  is  proportional  to  h"5  as  in  the  improved  Euler 

method.  The  three  term  Taylor  series  method  requires  the  calculation  of 

partial  derivatives  f  ,  f  .  In  principle,  even  higher  term  series  can 

x  y 

be  used,  but  the  difficulty  in  computing  the  terms  with  the  higher 
derivatives  makes  the  method  very  awkward  to  use. 

c.  Runge  Kutta  Method 

The  Runge  Kutta  method  involves  a  weighted  average  of  values  of 
f(x,  y)  taken  at  different  points  in  the  interval  xn  £  x  <_  xn  +  1 . 

A  commonly  used  Runge  Kutta  technique  is  characterized  by 

yn+l  =  £  [knl  +  2kn2  +  Zkn3  +  krV  (31) 

where 

k"i  ■  f<v  V 

kn2  «  f (x  +  yn  *  |  knj) 
kn3  -  f(x„  *  y„  ♦  |  kn2) 
kn4  =  f(xn  +  h«  yn  +  h  kn3) 


c 

The  local  formula  error  is  proportional  to  h  ,  and  thus  a  more  accurate 
formula  Is  obtained  at  the  cost  of  more  computation.  The  above  formula 
is  equivalent  to  a  five  term  Taylor  series  approximation;  the  "k"  factors 
substitute  for  the  need  to  compute  partial  derivatives.  This  method  is 
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one  of  the  most  widely  used  methods  for  the  numerical  integration  of 
ordinary  differential  equations  in  engineering  analysis,  and  is  one  given 
attention  with  the  collocation  method  in  this  project. 

3.  SUMMARY  AND  CHOICE  OF  METHODS 

The  Runge  Kutta  method  is  selected  over  the  Euler  and  Taylor  methods 
because  of  the  improved  formula  error  and  ease  of  use  over  series  cal¬ 
culations.  The  weighted  residual  methods  are  chosen  because  they  provide 
a  means  to  approximate  the  velocities  in  the  Euler  equations  by  a  polynomial, 
trial  function  process.  This  will  be  explained  in  more  detail  in 
Section  IV. 

4.  GRAMMEL'S  WORK  AND  MULTI  RIGID  ANALYSIS 

This  subsection  briefly  reviews  the  approximation  techniques 
attempted  by  Gramnel  and  the  computer  methods  known  on  multiple  rigid  body 
analysis.  These  methods  have  been  identified  in  the  literature  search 
stage  of  the  project.  They  are  not  directly  relevant  to  the  problem  under 
study,  but  their  prevalence  warrants  the  following  summary. 

a .  Grammel ' s  Work 

The  work  of  Richard  Grammel  (References  22-24)  is  summarized  here 
because  it  represents  a  twentieth  century  attempt  at  the  solution  of  the 
unsymmetric  top  problem  by  means  of  an  iterative,  trial  function  process. 


Grammel's  work  is  of  a  more  analytic  than  computation  nature. 
Grammel  considers  Equations  5-7  and  studies  trial  solutions  of  the  form 


t  Qt 

<*>2  =  a  + 

u>2  =  b  +  £2^^ 

<*>3  =  c  +  E3e^t 


(32) 


where  Q  is  a  complex  variable. 

Grammel  is  able  to  establish  stability  criteria  and  to  in  some  cases 
solve  the  problem.  The  approach  does  not  make  use  of  the  weighted 
residual  method,  and  the  mathematical  expressions  in  the  higher  iterations 
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become  cumbersome.  However,  foundations  are  laid  prior  to  the  introduction 
of  large  scale  digital  computers. 

b.  Multi-Rigid  Body  Analysis 

Multiple  Rigid  Body  analysis  is  an  efficient,  computer-oriented 
program  to  study  and  resolve  the  motion  of  complex  systems  of  linked 
rigid  bodies.  As  such,  it  is  not  relevant  to  the  rigid  body  problem 
under  study  here.  It  would  be  applicable  to,  for  example,  a  satellite 
with  flexible  appendages,  an  orbiting  antenna  composed  of  connected 
rigid  links,  or  the  bones  of  the  human  spine,  rib  cage,  neck,  etc.  A 
typical  chain  system  is  shown  in  Figure  2.  Finite  segment  and  multi- 
rigid  body  modeling  programs  contain  an  efficient  means  to  label  each 
link  and  its  connection,  then  to  reduce  the  system  to  its  equivalent 
motion  using  algebraic  operations  and  D'Alembert's  form  of  Lagrange's 
equations.  These  methods  are  in  essence  a  computerized  application  of 
planar  mechanism  theory.  As  such  they  do  not  address  the  issues  of 
interest  in  the  rotation  of  unsymmetric  tops.  However,  two  useful 
concepts  are  used  in  these  formulations.  One  is  the  use  of  hybrid  or 
quasi -coordinates,  which  is  discussed  in  Section  II  under  Lagrangian 
dynamics.  Second  is  the  use  of  Euler  parameters  or  quaternions  to 
eliminate  singularities  in  the  model. 
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SECTION  IV 

COMPUTER  SIMULATION  OF  THE  RIGID  BODY 

This  chapter  describes  the  application  of  the  collocation  method 
to  the  Euler  Lagrange  equations  (13-15)  In  a  trial  function  scheme. 

In  subsection  1,  the  rationale  for  trial  functions  Is  explained.  In 
subsection  2,  the  algorithm  for  use  of  collocation  Is  detailed.  In  sub¬ 
section  3,  the  algorithm  for  a  four  term  Runge  Kutta  method  Is  outlined. 

The  capability  to  model  problems  with  time  varying  torques  Is  added 
to  both  the  Runge  Kutta  and  collocation  programs.  The  Idea  to  include  this 
capability  results  from  a  survey  of  Grammel's  papers.  Examples  of  such 
torques  acting  on  a  top  result  from  a  rotating  electromagnetic  field,  a 
change  In  forces  acting  on  a  satellite,  or  the  influence  of  a  small  retro 
jet  on  a  space  body  when  the  percent  change  of  mass  is  Insignificant. 

1 .  TRIAL  FUNCTION  APPROACH 

The  Euler  Lagrange  equations  (13-15)  can  be  sunmarlzed  in  the 
following  form 

“1  •  Vj“k  =  Mi  (33) 

where  i,  j,  and  k  range  from  1  to  3  in  cyclic  order.  The  nonlinearity 

in  these  equations  is  caused  by  the  coupling  between  the  w.  u>.  terms. 

J  ” 

If  the  rigid  body  were  symmetric,  e.g.,  sphere,  all  such  terms  would  be 
zero  and  the  problem  would  reduce  to  the  simple  integration  of  =  M^ 
with  solution  =  M^At.  A  nearly  symmetric  top  would  have  coefficients 
A.j  nearly  zero  and  the  equations  might  be  modeled,  by  first  considering 
the  uncoupled  solution 

0>i  =  M.At 

to  obtain  an  Initial  estimate  of  the  terms,  then  re-adjusting  them 
to  fit  the  equations  with  the  coupled  terms  uj  formed  from  the  initial 
estimates.  Thus  a  linear  estimate  M^ At  could  be  used  for  each  Interval. 

The  problem  with  this  approach  Is  the  limitation  on  a  velocity  term  u^. 
which  Is  In  reality  changing  throughout  the  Interval  At.  In  other  words. 
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the  linear  estimate  M^At  allows  only  constants  for  u>^.  A  more  general 
top  would  require  that  the  intervals  At  be  made  small  enough  that  the 
constant  terms  still  track  the  solution.  In  practice  the  interval 
will  reach  a  lower  limit  and  the  approximation  will  fail  to  be  effective. 

The  choice  of  trial  functions  is  based  on  the  level  of  dynamical 
behavior  expected  in  the  problem,  and  by  the  computational  or  algorithmic 
constraints  on  the  programmer.  A  linear  top  problem  would  only  require 
linear  trial  functions  of  the  form  c^t  (c-|  =  M-j).  A  general  unsymmetric 
problem  would  require  higher  order  trial  functions.  In  this  effort  the 
choice  is  made  of  two  and  four  term  trial  functions  of  the  polynomial 
family.  The  polynomial  allows  the  derivative  terms  to  be  easily 
expressed.  A  four  term  trial  function  allows  i!>..  to  be  easily  expressed. 

A  four  term  trial  function  allows  <S.,  an  acceleration  term,  to  vary  up  to 

Cj  +  2c2t  +  3c3t2  +  4c4t3 

and  is  thus  a  rather  general  motion.  The  cost  in  programming  is  the 
expansion  from  the  solution  of  two  equations  in  two  unknown  parameters 
in  each  interval  to  the  solution  of  four  equations  in  four  unknown 
parameters.  Higher  order  trial  functions  can  be  written,  but  the  im¬ 
provements  are  limited  as  the  higher  powers  of  t  add  only  small  con¬ 
tributions  to  the  estimation. 

2.  COLLOCATION  PROGRAM 

• 

If  the  above  is  generalized  to  allow  at  least  a  first  order  change 
in  <*>..  throughout  the  interval  then  the  idea  of  trial  functions  follows. 
Consider  an  initial  estimate  of  given  by 

2 

=  0)^(0)  +  Cjt  +  c2t  (34) 

The  w^(0)  satisfies  the  initial  condition  on  u>. .  The  c-j ,  and  c2  are 
parameters  to  be  determined  in  the  interval  At.  Proceeding  with  an 
initial  estimate  of  ,  u)2,  o>3  the  u>j  terms  are  calculated  and  used  to 
solve  the  equations  for  the  interval.  The  reason  for  separation  of  the 
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problem  Into  a  first  estimate  and  a  correction  is  the  non-linear  coupling 
of  the  Euler  equations,  and  the  fact  that  u>j  are  not  available  to 
calculate  etc. 

To  illustrate  the  above  procedure,  let 
WjU)  =  u)j{0)  +  Cjt  +  c2t2 

u>2(l)  =  u>2(0)  +  d^t  +  dgt2  (34a) 

w-jd)  =  (o3(0)  +  ejt  +  e2t2 

This  allows  <^  =  c^  +  2c2t  =  a  linear  function  of  time  over  the  Interval 
under  study.  Now  form 

J>i  =  F  =  Mj  +  Aw2<i>3 

where  w2  w3  are  the  initial  estimate  and  are  found  from  the  initial 
conditions  a>2(0)  and  w3(0).  Upon  substitution  of  the  trial  functions  for 
di-| ,  w2,  i3,  the  following  expressions  result 

Cj  +  2c2t  =  F  =  Mj  +  A<o2(0)o»3(0) 

dx  +  2d2t  =  G  =  M2  +  B(i>j(0)b>3(0))  (35) 

ej  +  2e2t  =  H  =  M3  +  Cw1(0)w2(0) 


For  the  initial  estimate  a>-j  ( 1 ) ,  co2(l),  ui3(1 )  one  can  solve  each  of  the 
above  three  equations ^separately  by  expanding  each  by  one  of  the  four 
weighted  residual  methods  discussed  in  subsection  1  of  Section  III. 
Using  collocation,  two  points  within  the  interval  are  selected  and  each 
of  the  three  equations  is  evaluated  at  the  two  points  to  determine  the 
coefficients  c^ ,  c2,  ...  e2.  See  Figure  3  for  a  typical  interval. 

For  example 


Cj  +  2c2ti  =  F(tj) 

^1  +  2c2^2  =  ^(^2^ 

tj,  t2e(tQ,  t3) 


(36) 
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To  T-j  T2  T3 

I - 1 - 1 - 1 

Figure  3.  Time  Increment  for  Collocation 

For  the  Initial  estimate  It  Is  assumed  that  F  is  constant  over  the 
interval,  since  only  the  initial  values  of  w2,  w3  are  available.  The 
above  can  be  rewritten  In  matrix  form  as  follows 


1 

2tl 

V 

’F1 

1 

h 

2t2 

4 

C2 

F2 

k  , 

(36a) 


Having  determined  c-j  ...  e2,  initial  estimates  of  w1 ,  u>2,  w3  are  con¬ 
structed  from  the  constitutive  definitions  (Equation  34a).  These 

estimates  are  used  to  now  calculate  the  u>.  w.  terms  in  F,  G,  H  and  a 

J  K 

corrected  calculation  can  be  made  which  allows  these  terms  to  vary  with 
time  over  the  interval.  When  the  corrected  u.  has  been  calculated,  the 
process  iterates  forward  with  the  values  of  c^O),  o>2(  1 ) ,  u>3(l)  as  the 
new  initial  conditions. 

The  subdomain,  Galerkin,  and  least  squares  criteria  all  require 
that  an  integration  be  performed  on  each  interval  of  time,  whereas  the 
collocation  method  makes  use  of  the  sifting  property  (Equation  27)  to 
evaluate  the  residual  equations  at  specific  points  in  each  interval. 

The  integration  requires  that  the  terms  of  the  u>.  u>t  be  written  out  and 
and  evaluated  over  each  interval.  For  a  two  term  model,  one  such 
product  looks  like 

(cjt  +  c2t2)  (djt  +  d2t2) 

In  the  collocation  method,  one  need  only  to  calculate  each  ^  ,  uj2  and  u>3 
term  at  the  points  t-| ,  t2  within  the  interval  (see  Figure  3)  and  to  apply 
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them  In  the  u>.  u.  terms,  I.e.,  long  expressions  of  trial  functions 

J  * 

multiplied  together  are  not  needed. 

The  subdomain,  Galerkln,  and  least  squares  methods  were  Implemented 
in  Fortran  programs,  but  were  later  dismissed  from  further  consideration 
because  they  failed  to  produce  meaningful  results  for  the  test  cases  where 
the  solution  is  known.  It  Is  expected  that  the  computer  error  generated 
from  the  long  expressions  of  u>.  terms  caused  this  error. 

The  program  flow  for  the  collocation  method  is  shown  In  Figure  4. 

The  locations  used  in  the  two  term  models  are  1/3  and  2/3  over  each 
Interval;  for  the  four  term  model,  the  locations  are  .2,  .4,  .6,  and  .8. 
The  specific  algorithm  for  one  interval  procedes  as  follows: 


a.  DT  =  .01 

b.  T1  =  1/3DT,  T2  =  2/3DT 

c.  Solve  for  c^,  c2 


1 

2V 

f  * 

C1 

*  ' 

Fi 

+  Au2(0)u»3(0) 

1 

k 

2t2 

C2 

t 

F2 

k 

Mj  +  Ao)2(0)(.)3(0) 

and  likewise  In  dj,  d?,  e^  e2. 

d.  Form 

0^(1)  =  wjCO)  +  Cjt  +  c2t2 
»2(1)  =  w2(0)  +  djt  +  d2t2 
u>3(l)  =  «3(0)  +  ejt  +  e2t2 

e.  Now  calculate 

F1  =  =  Mi  +  Ao)2(t1)u3(t1) 

^2  =  ^(tg)  *  +  Aco2(t2)w3(t2) 


and  likewise  In  remaining  variables. 
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Figure  4.  Collocation  Program  Flow 
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f.  Solve 


* 

1 

2tl‘ 

r  ' 

C1 

F(tx) 

1 

2t2 

t 

c2 

b  4 

F(t2) 

*  * 

g.  Assemble 

Wj(l)  =  Ul(0)  +  c^t  + 

102(1)  =  002(2)  +  djt  +  djt2 
oo3(l)  =  (03(0)  +  e^  +  e2t4 

h.  Print  results  for  interval 

i.  Exchange  00. (1),  00. (0)  and  iterate  to  next  interval. 


3.  RUNGE  KUTTA  PROGRAM 

The  Runge  Kutta  method  from  Section  III,  Equation  31  is  implemented 
in  a  Fortran  program  for  comparison  with  collocation.  The  use  of 
Equation  31  is  straightforward;  the  Runge  Kutta  "k"  terms  are  calculated 
in  each  interval  and  the  estimates  for  00^,  ug*  003  are  assembled  from 
these  terms.  The  interval  size  for  both  the  Runge  Kutta  and  collocation 
is  chosen  to  be  .01.  In  comparing  the  Runge  Kutta  and  collocation  methods, 
the  same  step  size  is  used  throughout.  A  copy  of  the  Runge  Kutta  program 
can  be  found  in  the  appendix.  The  program  flow  is  illustrated  in 
Figure  5. 

The  Lipschitz  condition  is  used  in  the  Runge  Kutta  program  in  order 
to  guard  against  unbounded  approximations.  An  initial  Lipschitz  constant 
of  500  is  used;  in  some  cases,  larger  values  are  used. 
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Figure  5.  Runge  Kutta  Program  Flow 
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SECTION  V 
RESULTS 

The  collocation  programs,  as  well  as  the  Runge  Kutta  program  from 
Section  IV,  are  evaluated  for  effectiveness  by  using  a  number  of  test 
cases.  The  test  cases  represent  rigid  bodies  of  given  moments  of  inertia 
and  torques.  Subsection  1  of  Section  V  details  the  test  cases  used  in 
this  work.  In  subsection  2,  a  discussion  of  the  graphical  results  is 
included,  as  well  as  a  brief  discussion  of  error.  In  subsection  3,  a 
brief  discussion  is  provided  on  computer  execution  time,  complexity,  and 
general  concerns. 

1.  TEST  CASES 

The  following  test  cases  are  used  in  evaluation  of  computer  approaches 
from  Section  IV.  (See  Table  1.) 

Case  1  is  a  simple,  linear,  uncoupled  top,  and  is  included  for 
validation  of  the  methods.  The  solutions  are  constant  functions  of  time. 

Case  2  is  a  force  free,  symmetric  top,  and  can  be  found  in  Goldstein 
(Reference  1,  pages  161-162).  The  solution  is  known  and  is  given  by 

wj  =  sin  5t 

U)2  *  cos  5t 

w3  =  10 


Case  3  is  a  small  perturbation  of  the  problem  in  case  2. 

Case  4  is  an  unsymmetric  top  with  no  torques. 

Case  5  is  the  same  top  with  unequal  torques  added. 

In  Case  6  a  grossly  unsymmetric  top  is  modeled. 

In  Case  7  the  top  of  Cases  4  and  5  is  now  given  time  varying  torques. 

The  values  for  the  moments  of  inertia  and  torques  are  chosen  to  produce 
reasonable  numbers  on  output  and  have  no  special  physical  significance. 

In  Cases  1  and  2  an  objective  study  of  the  error  can  be  made  against 
known  solutions.  In  the  remaining  cases,  the  solutions  are  not  known  and 
one  can  only  compare  the  results  qualitatively. 
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TABLE  1 
TEST  CASES 


1. 

!1 

=  !2 

=  *3 

=  40 

N1 

=  n2 

=  n2 

-  100 

2. 

!1 

=  !2 

o 

II 

!3 

N1 

=  N2 

=  N3 

=  0 

3. 

h 

=  41, 

•  l2  = 

39, 

N1 

=  n2 

=  N3 

=  0 

4.  Ij  =  40,  I2  =  30,  I3  =  20 

N1  -  N2  -  N3  -  0 

5.  Ij  =  40,  I2  =  30,  I3  =  20 

Nj  =  100,  N2  =  200,  N3  =  150 

6.  Ij  =  200,  I2  =  10,  I3  =  20 
N1  =  100,  N2  =  200,  N3  =  150 

7.  Ij  =  40,  I2  =  30,  I3  =  20 

Nj  =  100,  N2  =  200,  N3  =  150 
(initially) 

Nj  =  sin  lOOt,  N2  =  cos  lOOt,  N3 
(as  t  >  0) 
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2.  GRAPHICAL  RESULTS 

The  results  of  the  seven  test  cases  are  plotted  in  Figures  6  through 
12.  The  following  is  a  discussion  of  these  results. 

In  Case  1,  Figure  6,  all  three  programs  (2  and  4  term  collocation, 
Runge  Kutta)  show  good  correlation  and  coincide  with  the  known  solution. 
There  was  no  deviation  in  the  results  of  any  of  the  programs  from  the 
known  solution. 

In  Case  2,  Figure  7,  the  known  solution  is  plotted  along  with  the 
2  and  4  term  collocation  and  Runge  Kutta  programs.  The  collocation 
programs  show  a  good  correlation  with  the  known  solution  and  are 
coincident  with  the  known  solution  for  the  graph  shown.  The  error  was 
found  to  be  of  the  order  of  .05  to  1%  and  constant.  The  Runge  Kutta 
program  gave  poor  results,  and  terminated  after  only  25  steps  due  to 
the  Lipschitz  condition.  The  w-j  curve  can  be  seen  to  be  growing 
beyond  the  expected  sine  solution.  The  Runge  Kutta  program  has  been 
rechecked  for  an  error,  but  none  has  been  found. 

In  Case  3,  Figure  8,  the  solution  of  Case  2  is  slightly  perturbed 
such  that  w3  will  no  longer  be  a  constant.  The  collocation  programs 
concur  with  one  another.  The  Runge  Kutta  program  is  again  diverging. 

In  Case  4,  Figure  9,  an  unsymmetric  top  is  modeled.  The  collocation 
programs  concur,  but  the  Runge  Kutta  program  fails. 

In  Case  5,  Figure  10,  one  notices  some  differences  between  the  2  and 
4  term  collocation  models.  The  4  term  model  is  assumed  to  be  the  more 
reliable  result  because  of  the  extra  terms  and  the  experiment  described 
in  Case  6.  Figure  10  is  divided  into  an  A  and  B  part  to  illustrate  the 
different  curves. 

In  Case  6,  Figure  11,  one  sees  major  differences  between  the  models. 
The  2  term  model  diverged  and  is  not  shown.  The  4  term  model  illustrates 
very  dynamic  results.  Because  of  the  disparity  between  the  2  and  4  term 
results,  it  is  judged  that  the  motion  being  simulated  is  too  fast  for  the 
2  term  model.  To  check  the  4  term  collocation  model,  3  and  5  term  models 
have  been  used.  The  results  of  the  3,  4,  and  5  term  collocation  models 
concur. 
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In  Case  7,  Figure  12,  the  2  and  4  term  models  concurred  fairly 

well . 

3.  SUMMARY  COMMENTS 

The  collocation  programs  use  approximately  twice  the  execution 
times  of  the  Runge  Kutta  program  to  compute  100  steps.  However,  this 
is  of  no  special  concern  as  the  average  execution  times  are  of  the 
order  of  .3  seconds,  which  is  negligible  by  current  computing  standards. 

The  collocation  programs  are  of  about  the  same  level  of  complexity 
to  program  and  use  as  the  Runge  Kutta  method.  The  collocation  methods 
require  an  understanding  of  the  use  of  a  matrix  subroutine. 

The  step  size  has  been  changed  from  .01  to  .001  for  several  cases 
in  order  to  evaluate  error  trends  in  the  programs.  The  collocation 
programs  provide  generally  better  results  with  this  step  size  (with 
the  exception  of  Case  6,  2  term)  whereas  the  Runge  Kutta  model  does 
not  improve. 
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2  TERM  COLLOCATION 
4  TERM  COLLOCATION 
RUNGE  KUTTA 


TIME  SECONDS 


Figure  6.  Case  1 


KNOWN  SOLUTION 
2  TERM  COLLOCATION 
4  TERM  COLLOCATION 
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Figure  7.  Case  2 
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.  2  TERM  COLLOCATION 

-  4  TERM  COLLOCATION. 

- RUNGE  KUTTA 


.6  .7  .8  .9  1.0 


TIME  SECONDS 


Figure  8.  Case  3 
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Figure  9.  Case  4 
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0  .1  .2  .3  .4  .5  .6  .7  .8  .9  1.0 


TIME  SECONDS 


Figure  10.  Case  5 
PART  A 
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Figure  10.  (Concluded) 
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TIME  SECONDS 
Figure  11 .  Case  6 
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2  TERM  COLLOCATION 
4  TERM  COLLOCATION 


Figure  12.  Case  7 
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SECTION  VI 

SUMMARY,  CONCLUSIONS  AND  RECOMMENDATIONS 

The  general  rigid  body  problem  of  subsection  1  of  Section  I  has  been 
modeled  by  a  trial  function  technique  on  a  digital  computer.  The  trial 
functions  represent  estimates  of  angular  velocities  of  the  rigid  body 
over  intervals  of  time.  The  results  of  the  computer  models  have  been 
studied  against  known  solutions  and  graphically  against  another  numerical 
method. 

The  results  of  Section  V  allow  a  number  of  conclusions  to  be  made. 

The  weighted  residual  method  is  a  successful  trial  function  approach  to 
the  rigid  body  problem.  In  particular,  collocation  is  the  successful 
technique  of  the  weighted  residual  methods.  Collocation  is  superior  to 
the  Runge-Kutta  method  employed  In  this  project.  Collocation  is  of  about 
the  same  level  of  computational  complexity  to  program  as  Runge-Kutta. 

The  results  indicate  that  a  three  or  four  term  collocation  model  is 
probably  sufficient  for  most  problems. 

The  above  conclusions  leave  a  number  of  issues  open  for  further 
Investigation.  The  following  recommendations  are  suggested  for  further 
work: 

1.  Experiment  with  alternative  step  sizes  and  trial  functions  for 
collocation,  and  study  practical  matters  of  convergence  and  error. 

2.  Extend  the  collocation  method  to  a  predlct-correct  technique  and 
evaluate  against  prior  results. 

3.  Compare  the  results  of  this  project  with  another  well  established 
numerical  Integration  routine,  e.g.,  Adams  predictor  corrector. 

4.  Use  the  results  from  above  to  Investigate  classical  solutions  in 
the  literature. 

5.  Develop  a  means  for  transforming  ,  u>2*  “>3  from  body  coordinates 
to  space  coordinates  and  Integrate  to  Euler  angles. 
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APPENDIX  AND  PROGRAMS 

2  Term  Collocation 
4  Term  Collocation 
Runge  Kutta 


AFWAL-TR-81-1245 


L,A 


IN* 

m* 

1 29* 

IN' 

IN' 

IN' 

1S4-C 

174-C 

1M'C 

194*C 

set* 

214* 

224- 

234'C 

244*C 

254* 

244* 

2?4« 

2t4'C 

294  *C 

344* 

314* 

324- 

334* 

344*1 

344* 

344*4 

374* 

344 *C 

394*C 

444* 

414* 

424' 

434*C 

444*C 

4S4* 


474* 

444*C 

494-C 

544* 

514* 

424* 

S34*C 

S44*C 

544* 

S44- 

474* 


PROGRAM  FEN1  <  INPUT , OUTPUT ,  TAPES >  TAPES ) 

DIMENSION  4(2, 2). I <2, 4), US* (12) 

REAL  I1,I4.I3»N1.N2,N3,M1,N2.N3.NS1,NS2.NS3 
REAL  UP 

DIMENSION  Ul(4),g2(6),M3(4),Ul(4).W2(S',U3C«) 

TINE-4. 

COLLOCATION  PROQAAM 
TWO  TERN  APPROXIMATION 
HOMCNTS  OF  INERTIA 
11*44. 

12*44. 

13*24. 

TOROUE4 

Nl*4. 

N2*4. 

N3*4. 

INITIAL  CONDITIONS  U 
(ltd  1*4. 

U2(t )*1. 

U3(l)*14. 

URITE(6f  1  )I1,N1,I2,N2,I3,N3 

F0«NAT(2X.FC.2.SX,FS.2,/,2X,FS.2,SX,F6.2./,2X,FS.2.SX,F6.2) 

URITK4.4) 

FORMAT* 6X,2HU1, 14X.2HU2, 1RX,8HU3, 16X, IHI, 16X.2HU1, 
C24X,2HU2,24X,2HU3) 

NORMALIZED  MOMENTS  OF  INERTIA 
A1*(12-I3)/I1 
»l*(I3-Il)/I2 
C*<I1-I2)/I3 

NORNALIZED  TORQUES 
N1*N1/I1 

N2-N2/I2 

M3-N3/I3 

INITIAL  VELOCITIES 
Ul(l 1*4. 

U8(l )*4. 

U3<1)*4. 

MASSES 
MSI *144. 

M48*114. 

MS3*184. 
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sm*c 

CM*C 

FORCES 

(M> 

FRC1-1C. 

C1C* 

F1K2-H. 

Ut- 

FRC3-4C. 

CM* 

ACL1*FRC1/NS1 

MI¬ 

ACL3-FRC2/NS2 

CH- 

ACL3*FRC3'NS3 

CM* 

LIF*SM. 

c?c* 

DT..C1 

CM*C 

CMC  INITIALIZATION 

CM* 

M  IN  1*1,333,1 

7M* 

TIH£*TINE*DT 

710* 

Tl-lT/2. 

7M- 

T2-e.*ftT/3. 

7M* 

TS-.Ci 

74C*C 

RIGHT  MANS  SIDES 

7M* 

F1*N1*AI*U2(1)*U3<1> 

7M* 

Gi*n24»inini)tu3(n 

770* 

M1*N3*C3U1 ( 1 )SU2< 1 ) 

7M* 

Xl-ACU+uai  1  )*U3<  1  )-V3(  1  )SU2(  1 ) 

7M* 

VI • ACL2*V3  ( 1  J3U1  ( n-Ul  ( 1  >SU3  ( 1 ) 

CM* 

Zl • ACL3*U1  ( 1  )*U2(  1  )-U2(  1  )SU1  ( 1 ) 

•1C*C 

COEFFICIENT  MATRIX  A 

C2C- 

A(l,t)*l. 

CM* 

A(1,2)*2.»T1 

Ml* 

A<2.1)*1. 

CSC* 

At  2,2 )*2.XT2 

CM*C 

FORCING  FUNCTION  MATRICES 

C7C* 

B(l.l)*Fl 

CM* 

3(2,1 }*F1 

CM* 

3(1.21*01 

MC* 

3(2.2)*G1 

91B* 

K1.3)*H1 

B2C* 

B(2.3)*H1 

CM* 

C(1,4)*X1 

MC* 

»<2.4>*X1 

CSC* 

»<1,S)*V1 

CM* 

»<«,S)*V1 

C7C* 

3(1.0*21 

CM* 

9M*C 

3(2.0*21 

MATRIX  SUBROUTINE  PARAMETERS 

1CCC-3C 

CONTINUE 

ictc* 

N*C 

icac* 

N*2 

IBM* 

IA-2 

IC4C- 

IDOT*S 

1CM* 

CALL  LEGTlF(A.n.N.IA.3,IDCT.USAl 
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154** 

tSM* 

1W> 

1S?M* 

ISM* 

l»li* 

!!«• 

ISM* 


ISM* 

ISM* 


1M 


ut(4)-ui(i  H»a,5  >«♦•<«.«  )rrsm 
v3(4)>u3(n«s(i,s)m«i(s,s>tTtttt 
MSm<S.S>Ul(4),y«(4),U3(4)tI,ICT,UlC4),US<4>,W3(4) 
f0WMr(ix,ris.4.iSK,ris.4.iM,ris.4.SK.ia.».ii>M#Fis.4. 

C10X,F10.41HX,F1#.4) 

UftlTt(S>TIHC,Vl(4),U8<4),U3(4  )>UK4>,UB(4),U3(4) 

UKl  >*W1<4) 

VIU)<Ut(4) 

U3(l )-U3(4) 

V1(1)*U1(4) 

U2(1)*UR(4) 

ui(n-u3(4) 

OT-0T+.S1 

CONTINUE 

END 
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U*lO* 

\X 

13* 

14* 

IS* 

1C* 

m«c 

1S*C 

IQQ-C 

acc-c 

ci* 


M*C 

«*C 


r». 


IM-C 

3QQ-C 

310* 

SS' 

35#  ■  1 

3M« 

37*4 

31* 

39*C 

4S*C 

41* 

48Q* 

43*. 

440*C 

4S*C 


47C- 

4M- 

4QQ-C 

CM-C 

11* 

SC* 

S3* 

S4*C 

ssc-c 

sc* 

S7Q. 


mum  «i»iai#oT.ouTcut,rij|s,w5ui 

dihehsion  *<4,4>,»c4,c»,usiM«4>,T<i>,r<i>.c<s>,m»>. 
ccMt<ii|ii«ra!Ni.NC,N3<m.M.ia.Mi.Ma.M3 
M&NSION  Ml<«Mtf<l)#U3tSl,Ul(Q>.Ut<C>.U3(f> 

nii*#. 

COLLOCATION  MOOMN 
roue  TUN  ATMOXlRATION 
NOHENTS  or  INERTIA 
I1-4C. 

IC-4C. 

I3-M. 

TORQUES 

m*Q. 

NC«Q. 

N3«C. 

INITIAL  CONDITIONS  U 
U1  <!>•#. 

U2(l)*l. 

U3U  >»1S. 

roR«AT(8x.|Fs!S!s5ff*a!^S§»r*.c.SN,F*.c,/’,ax,r«.c»sx,rc.c> 

FORNArfcx/aHui ,  lsx.axua,  itat.ma.  wx,  in:  ,  isx.aHui . 
caox.EHua.asx.aHws) 

NORHALIZED  NONCNTS  OF  INERTIA 
Afda-I3)/Il 
si*(i3-in/ia 
c*di-ia>/i3 

NORfMLIZED  TORQUES 

m-Ni/ii 

H2-N2/I2 

H3-N3/I3 

INITIAL  VELOCITIES 
Ul(l>-0. 

ua<i)-c. 

U3<1 >■•. 

HASSES 

HS1-10Q. 

RS8*110. 


AFVJAL-TR-81-1245 


■if 


OS#- 


•1#*C 

S: 


men 

acu*fnci/hsi 

ACU*FRCI/NSt 

KU'FKl/m 

L1R-S00. 

OT-.it 

CN»  INITIALIZATION 
00  1M  1*1,300, 1 

riMc-riMEor 

n-oro.t 

trots.  4 

T3-0T0.S 

T4OT0.8 

T5-.01 

aiomt  mm  sioes 

Fl«Mt4AltU2(l )SU3(1 > 

ci«H2oiavm>sii3(i> 

H1»N3»C*U1(1  XUSCt  ) 

XI •ACll*V2( 1 )0V3< 1 )-03( 1 )SU>( 1 > 

Y1-ACL2+V3U  mmi  Hill  )«I3(1 ) 
zi*ACL3+wiu)a«ci)-vaci)«ii<i  > 
COEFFICIENT  matrix  a 
A(l,|  )■}, 

A(l,2)*e.*Tl 

AU.3)*3.«T»*2 

A<1,4)*4.*T1**3 

A<a,n*i. 

A(2,2)*2.*T2 

A<2,3).3.«T2»*2 

A(2,4)*4.*T2«3 

A<3.1)-1. 

A(3.2)-2.*T3 
A(3,3).3.*T3W2 
AO,4  >»4.*T3**3 
A(4,l)*l. 

A(4,2)»2.tT4 
AM,3)»3.3T4*S2 
A<4,4 J«4.ST4**3 
FORCINO  FUNCTION  MATRICES 
S(1,1)*F1 
•<2,1)01 
•0.1)01 
0C4,1).F1 
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00  40  I't.S.t 
Ttaj-ri 
TO  NT* 

TC4NT3 

T<SNT4 

wt  a  nui  ( i  )«•(  i .  i  >*T(i  >♦»<*,  Dtrct  >«**♦*<  3.  t)*ra  >«»♦ 

0<4,l)*TU)t(4 

M(L)«us(iH>u.a)>Tcu«»(*#*»ra>m4»o.*»Ta»04 

ci(4(i)ira)H4 

W3<lNV3UN»<l,3>*TaN»<*,3>*T<l>MM(3.3>rr<l)M* 

0(4,3J*T(l)W4 

VI  (L  NVK1  )♦•(  X  ,4  )*T(l  NKI.4  )tTO.  >tt*4»(3,4)tm  )U3* 
0(4.4)*m»>4 

u*a  )«v*(  i  )♦»<  i  ,b  >mu4t<*.*  >rra  ittMtu  wm  »*»♦ 

0(4.»  irriuw 

W(tNV3ciN»(i,»>rra>4»<*,*)tT(um4i<3.«)*Ta>M3* 

0<4.«>ST(l)M4 

F(L)«Rl-HUtU*(U<U3(L) 

•a  nm+si«u  a  ma  ci > 

H(u«R3*cnit  a  )«uta ) 


AFWAL-TR-81-1245 


>SM« 

lf?2M22 

ism-c 


x<u*«eu«uitL  mint  >-v3(i>*u«<u 
V<  L  )>AC12*U3  (l)ttll(l  >-Ul <1 >*W3( l ) 
Z(  L )  *<*Ct34WI  (t )«!*( IM*<  L  )«W1  ( U 
CONTINUE 

NOW  KOO  INTEftlMl 

au,d«i. 

A(1,2)*«.ST1 

A(l,3)»3.*Titl2 

*(t,4)<4.tun3 

A(2,2)*2.ST2 

A(2,3)«3.*T2M2 


A(2,4)>4.«T2ST3 

A(3,l)-1. 

A(3,2)*2.IT3 

A(3,3)*3.*T302 

A(3,4)*4.»T3»*3 

A(4, 1 )*1 . 

A(4,2)»2.*T4 

A(4.3)»3.*T4M2 

A(4,4)»4.*T4«*3 


»(1,1>-F(2> 


1722-  1(2,1  )*F(3) 
1722-  1(3, 1  )*F(4> 
IBM*  1(4, 1  )*F(5 ) 
1210-  1(1,2)*C<2) 
1222*  1(2.2  >*G( 3 ) 
1130*  1(3,2)*G(4) 
1240*  »(4.2>-G(5> 
1252*  1(1,3)-H(2) 
1222*  l(2,3)*H(3) 
1272-  1(3.3)-H(4) 
1222*  B(4.3)«H(5) 
1292*  1(1 , 4  )»X<2 ) 


2212* 


l(2,4)*X(3) 

»(3.4)-X(4J 

»(4.4)-X(5) 

•<1.5)*V(2> 

K2,5)*V(31 

»(3.5)-V(4) 

2(4.5)-V(5) 

>(l,OZ(2) 

l(2,S)*2(3) 

2(3.6)«Z(4) 

2(4, S )*Z(5 ) 

CAU  LEQT1F(A,N,N, IA.I, I0QT.U2A, ICR) 
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KAIL 


I##* 

no* 

ui> 

13#* 

14#* 

1M<C 

lf#*C 

17#* 

IM*G 

19#* 

I##* 

at#* 

2t#*C 

ta#-c 

84#* 

as#* 

as#* 

87#* 
■I 

_»*C 
3##*C 
31#* 
32#* 
33#* 

34#*C 

3S#*C 

38#* 

37#* 

38#* 

39#*C 

44#  *C 

41#* 


43#* 

44#*C 

48#*C 


47#* 

4M* 

4##*C 

6##*C 

ft#* 

ia#* 

S3#* 

54#*C 

554- 

f##- 

IT#- 


MAI  It«  18, 13«Nl,Nt.N3.AljN8.«N3,Nfi,Nll*N#J 
MAI  SF(4),IC#<4),KM<4>#ICX<4  i.KV(4),KZ«4» 

DINENsloN  UI<4MII<4MI3C4>,MU4>,U8l4>.M3t4) 

RUNQE  ICUTTA  F8088AN 

HOMNTf  Of  INERTIA 
11*4#. 

18*3#. 

13*8#. 

TOROUCf 

Hl*t##. 

N8-8##. 

URITE<f»ttIt*Nl,I8»N8«I3*N3  r<  I  tv  n  1  /  M  r|  l.fX.Ff.8) 

INITIAL  CONDITION#  U 
Ul( l  )*#. 

U2(l >*!• 

U3tl)*t. 

NORMALIZED  HOMNTf  OF  INMTIA 

A*(ia-i3vn 

D*ci3*nvia 

c*tii-ia»'i3 

NORMALIZED  torques 

M*N1/U 

N8-N8/I8 

N3-N3/I3 

INITIAL  WELOCITIE# 

UKI  )•#. 

U8(l )•#• 

U3(l )•#• 

MASSES 

ftfl*t##. 

M8*U«. 

N53M2#. 

FORCES 

FRC1*1#. 

FRC3-4#. 

FRC8*M. 
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m 


lilt* 
lilt* 

ii$ 

II4I- 

1IN* 

t# 

H7t* 

UN* 

lltt* 

IjT.: 

itM* 

UN* 

1240* 

12M- 

1M0*CC 

i*?t* 

lWt- 

13tt* 

131#* 

1310* 

1330* 

JM»- 

136*- 

136f- 

1371* 

1310* 

139*- 

14tt*C 

Hit* 

net* 

1430* 

144#* 

Hit* 

Htt* 

1470-C 

14«t* 

Hit* 

IMt- 

Iflt* 

ittt- 

ll3t- 


w3(i)*t3(j>4|ra(N(i> 

vuiwtUHVtmtf.iV 

mcn«wj>4imv(j> 

1*4 

J*3 

nr  1 1  >«m«MNR<  j  >*0(  j  y 

n(l)*«4MW3(J»WilJ) 

KN(I)*IQ4C<UiCJ)MK(J) 

0(1  I  )-tCU«Wt(  J  >«M3<  J  >-V3(  J  >«(  J  > 

KVC  i)*ACtl«V3(  J )  Hit  (  J  )-0i  ( J  >MO(  J  i 

n<t)'*ci3+uuj)nm(J)-**(J)**n(J> 

Ui(l)*«Mi«(OT/|.  )I(KF(1  >-*2.«(KF(B)**F (3) >*t(F(4> 
tft(I)-UM*(OT/|.  )S(K6(1  )*I.*(W(D'HCC(3 )  )*KQ(4 ) 
U3(I  )*UU3«(0T/C.  )f  (KH(1  Hl.KKMiOHCHO)  )«KH(4) 
01(1 >*001 ♦(0T/«.  )«CKX(1>*».*(W(2)*ICX(3))*CX(4) 
08(1  >*UU8*(0T/8.  )>(KV<1  H8.t(FY(2)*KV(3)  HKYM) 
w3<n*w3+<0T/t.  luau  >+a.t<a<8)*«Z(3  n+tz<4 ) 
URITt(t»S)Ul(I )»U2< I )<U3(Z)<L<Vi(I ),W81I )»W3< I ) 
IMlTCtSiTine.UKI  >.M2(1  ).W3CI  >,OKI  >.08(1  ).U3<I ) 
LIPSCMITZ  CONDITION! 

FL2*AI$(U1(I  MMl ) 

CU*A»5(U2(IMM8) 

HU*AIS(U3(I>-UU8) 

XL8*A»$<01(1)-001> 

VL8*A«f(W8(I>-008> 

ZL8*A»$(U3(I)-OU3) 

FU*AIS(ICF(4>-ICF(1>> 

GL1*AIS(KC(4)-KQ(1 ) > 

HU*ABS(KM<4>-KHU» 

XU*A«S(KX(4)HCX(1>) 

vu*aisckvm)-icv(ih 

ZU*A»S(KZ(4)-«Z(l)> 

LIP*Stt. 

IF(X11.GT.XL»IIP>00  TO  I 

IF(Yll .CT. YL23LIP >00  TO  I 
1F(ZL1.CT.Z12«IP)C0  TO  I 
IF(FL1.GT.F18«IIP>00  TO  I 
IF (CL  1.GT.GL8XLIP >00  TO  I 
IF(HU.QT.HL»LZP)60  TO  • 

UW1  *  U1(I> 

uua  •  u2( i ) 

IMS  •  13(1  > 

001*01(1) 

008-08(1 ) 

003*03(1) 


